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Abstract. We use three-dimensional numerical simulations to study self- 
organization in supersonic turbulence in molecular clouds. Our numerical exper- 
iments describe decaying and driven turbulent flows with an isothermal equation 
of state, sonic Mach numbers from 2 to 10, and various degrees of magnetiza- 
tion. We focus on properties of the velocity field and, specifically, on the level 
of its potential (dilatational) component as a function of turbulent Mach num- 
ber, magnetic field strength, and scale. We show how extreme choices of either 
purely solenoidal or purely potential forcing can reduce the extent of the iner- 
tial range in the context of periodic box models for molecular cloud turbulence. 
We suggest an optimized forcing to maximize the effective Reynolds number in 
numerical models. 



1. Introduction 

Modern statistical theories of fragmentation of molecular clouds (MCs) and star 
formation are based on an interpretation of the non-thermal emission linewidths 
and their correlation with length scale in terms of supersonic turbulence (Kaplan 
& Pronik 1953; Larson 1981; Heyer & Brunt 2004). It is believed that both 
the star formation rate and the initial mass function of newly born stars are 
controlled by MC turbulence (e.g., Padoan et al. 2007; McKee & Ostriker 
2007; Padoan & Nordlund 2009). There are competing views on the origin of 
this turbulence, which is either explained as a transient phenomenon associated 
with the cloud formation process or as if it were continuously driven by various 
energy sources (e.g., differential rotation, supernovae, stellar winds, protostellar 
outflows, variable FUV background, etc. (|Mac Low Klessenll2004l )). Since the 
typical Reynolds numbers in MCs are ~ 10 8 , expectations to find purely laminar 
regions in the cold star forming molecular gas on scales from ~ 50 pc d own to 
a few astronomical units should be pretty low (jElmegreen Sz Scald" 12004! ) . 

The supersonic regime typical of MC turbulence is extremely hard to achieve 
in the laborato ry and the informati on available from astronomical observations 
is limited (e.g-. lHever &: Bruntll2004i ). Most of what we know about the statistics 
of supersonic turbulence comes from large-scale numerical experiments intended 
to reproduce the basic non-linea r processes operating in the energy cascade in 



the inertial range of scales (e.g., iKritsuk et al.l l2007al ) . The effective Reynolds 
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numbers normally achieved in such simulations are at most ~ 10 4 , i.e. much 
smaller than the realistic values. Since computational resources are always lim- 
ited, even with the most advanced and least dissipative numerical methods only 
a short stretch of t he inertial interval can be captured at current grid resolutions 
up to 2048 3 zones ( Kritsuk et al.l l2009al ) . Simulations also rely on a number of 
assumptions intended to simplify the model and make it more tractable, such 
as the periodic boundary conditions and an isothermal equation of state. To 
achieve a better scale separation, a n approach of implicit large eddy simulations 
(ILES) is used ( Svtine et al.l 12000). The effects of molecular viscosity are, thus, 
replaced by numerical diffusivity of purely artificial nature. In simulations in- 
volving magnetic fields, the magnetic diffusivity is also replaced by the effective 
one built into an ideal MHD solver in use. In these circumstances the effec- 
tive magnetic Prand tl numbers usually achieved are on the order unity (e.g., 
Kritsuk et alll2009bh . 



The assumption of isothermality naturally restricts the physical box size 
to L < 5 pc, as the presence of multiple thermal phases plays a role on larger 
scales. In these circumstances, an artificial stirring force is required to mimic the 
turbulent energy flux from larger-scale cascade that cannot be modeled directly 
due to a finite physical dimension of the computational domain. With some 
rare exceptions, most of the models rely on a purely solenoidal forcing on large 
scales to provide a bett er 'boundary conditio n' for the inertial interval in the 



wavenumber space (e.g.. lBoldvrev et al.ll2002al ). While the stellar energy sources 



would mostly generate compressive fluctuations, the choice of a solenoidal force 
was then justified by the small compressional-to-solenoidal ratio measured in 
simulations. 

With higher quality and larger simulations available today, we can now re- 
assess the domain of applicability of solenoidal forcing in isothermal simulations 
of MC turbulence. To achieve this goal, we use various simulations we have per- 
formed in the past to investigate self-organization in supersonic turbulence and 
quantify the equilibrium compressional-to-solenoidal ratio in the inertial range 
as a function of the sonic and Alfvenic Mach numbers. We then come up with 
an optimized prescription for the large-scale forcing in isothermal periodic boxes 
that allows us to achieve a better scale separation. We show that a poor choice 
of forcing parameterization can easily lead to a complete elimination of the in- 
ertial range and result in rather "pathological" statistics, which have nothing to 
do with turbulence that fully develops far from boundaries and external forces 
at very high Reynolds numbers. 



2. Dilatational and solenoidal motions in supersonic turbulence 

The velocity field u(x, t) can be decomposed into solenoidal and dilatational 
parts u s and u c , such that u = u s + u c , V • u s =0 and V x u c = via Helmholtz 
decomposition. Let us consider = P(u- C ,k)/P(u s ,k) as a measure of the 
flow compressibility^ where -P(a, k) is the three-dimensional power spectrum of 
a vector field a, and k is the wavenumber. In an incompressible fluid x{k) = 0, 



1 There are also alternative measures of compressibility, e.g. 7 = (wj?) / (its), which gives the 
global ratio of the specific kinetic energy contained in the dilatational and solenoidal modes; 
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while in a compressible gaseous medium with purely potential (rotation-free) 
velocity field x(k) = oo. 

In fully developed supersonic isothermal turbulence, x(k) describes a bal- 
ance established via nonlinear exchange between the dilatational and solenoidal 
modes, which is ultimately controlled by the sonic (M s ) and Alfvenic (Ma) 
Mach numbers. In non-magnetized flows (Ma = oo) at low sonic Mach numbers 
(M s <C 1), compressibility is very weak and x(k) tends to settle at zero. At high 
turbulent Mach numbers (M s 3> 1), the compressional-to-solenoidal ratio hov- 
ers a round 1:2, which can be explained by simple geometrical considerations 
(e.g., iNordlund &: Padoanll2003l ). In MHD turbulence, the natural tendency to- 
wards Alfvenization, or dynamic alignment between the velocity field u and the 
divergence-free magnetic field B in the bulk of the volume away from shocks 
and dynamic rarefactions, results in suppression of dilatational activity. Thus, 
the presence of dynamically important magnetic fields would effectivel y reduce 
x(k) in trans- and sub- Alfvenic turbulence (e.g., Boldvrev et al.|[2002al ) . 

In simulations, besides M s and Ma, the ratio x(k) would also depend on 
the content of solenoidal and dilatational modes in the large-scale forcing, Xf- 
If X/ is fax from the equilibrium ratio that corresponds to chosen values of 
M s and Ma, the effect of such forcing will be felt further down the hierarchy 
of scales and the inertial range in such simulations would shrink or disappear 
depending on what Xf' ls enforced at the driving scale kt. Most of the simulations 
conservatively used a purely solenoi d al forcing with xt = , with the exception 
of Xf ~ 0.7 in lKritsuk et all (|2007aft . ISchmidt et "all (|2008l ) recently considered 
both Xf = and Xf = °o m 1024 3 non-magnetized simulations at M s ~ 5.5 and 
found that the velocit y scaling varies substantially with the large-scale forcing. 
Federrath et al.l (2008) also discovered that the density pd f in their compr e ssivel y 



driven models d o es no t bear a lognormal shape, see also ISchmidt et al.l (2009) 
iFederrath et~all (|2009h showed that x(k) ~ 1-2 at k/k min 6 [3,70] in their 
1024 3 simulation at Xf = oo, while at Xf = they obtained the expected 
X(k) « 0.5 at k/k min G [8,30], where k min = 2ttF1 ISchmidt] (120091 ) found that a 
transition from xj = 0tox/ = °o causes strong variations in spectral properties 
of turbulence in his large eddy simulations (LES). While these "pathological" 
statistics observed in simulations with purely compressive forcing clearly indicate 
a complete absence of an inertial range even at a grid resolution of 1024 3 zones at 
Xf = oo, they also hint at a possibility to optimize the problem setup by tuning 
the forcing to match the expected statistical equilibrium in the inertial range 
determined by the flow parameters. This would help to maximize the extent of 
the inertial range and thus to provide higher effective Reynolds numbers at the 
same computational cost. 



Xc(fc) = P(u c , k) / P(u, k), which estimates the fraction of dilatational modes in the velocity 
power spectrum as a function of k\ and r ca = (|V ■ u| 2 ) / ((|V ■ u| 2 ) + ^|V x u| 2 )), which 
represents the small-scale compressive ratio. Both \c and r cs are bounded in the interval [0, 1], 
while x an d 7 can potentially take arbitrary positive values. 

2 The same simulations also produced an unusual peak in x(k) on small scales at k/kmin G 
[300, 400] with the peak values of 2.3 and 4.0 in the runs with solenoidal and compressive 
forcing, respectively. 
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Figure 1. Statistics of supersonic turbulence from simulations with PPM 
and PPML: (a) velocity power spectrum, and spectra for dilatational and 
solenoidal parts, M s = 10 and Ma = 3; (b) x(k) for simulations with mixed 
forcing (xf = 0-7, red line) and with solenoidal forcing (\f = 0, all the rest); 
(c) pdfs of the alignment angle for three 512 3 simulations with Ma — 1, 3, 
and 10; (d) convergence of the cosf? pdfs for runs at 256 3 , 512 3 , and 1024 3 ; 
(e) same as (b), but for logiox(k) at 1024 3 only; (f) x(k) for 10 snapshots 
from a 1024 3 simulation of turbulence decay; (g) convergence of x(fc) for the 
5th snapshot from 512 3 and 1024 3 decay simulations; (h) x(k/k m i n = 1) as a 
function of M^ from a 1024 3 simulation of turbulence decay. 
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To explore the effects of artificial large-scale forcing on the velocity field 
statistics at scales adjacent to the forcing range, we collec ted data from various 
isothe rmal simulations with and without magnetic fields (jKritsuk et al.l [2007a . 



2009sjJ3) • The nonmagnetized runs u tilized the Piecewise Parabolic Method 
(PPM) of IColella fc Woodward! (11984T ) implemented in the ENZO code0 The 



MHD simulations wer e carried out with ou r Piecewise Parabolic Method on a 
Local Stencil (PPML. lUstvugov et al.ll2009l ). 

Figure la gives an example of the velocity power spectrum P(u, k) and 
products of Helmholtz decomposition for a 1024 3 simulation with Xf = 0j 
M s = 10 and M A = 3 (iKritsuk et al.1 l2009bh . One advantage of MHD sim- 



ulations, even in the super-Alfvenic regime, where the weak field is dynami- 
cally unimportant in most of the simulation domain, is the absence of visible 
bottleneck contamination in the inertial subrange adjacent to the dissipation 
range. This simplifies the discussion of the inertial range scaling for x(fe). In 
Fig. lb we collect the x(k) functions from various simulations. For instance, 
the red and green curves represent tw o non-magnetized sim ulations with PPM 
at M s = 6 with X f = 0-7 (1024 3 , iKritsuk et all l2007aD and with X f = 



(2048 3 , IKritsuk et al.1 [2009al ) . In both cases x{k) is close to the asymptotic 



value of 0.5 at wavenumbers k/k m i n G [10,100], as expected. Similar x-levels 
were also achieved in solenoidally driven non-magnet ized simulations at M s > 5 
by others ( Pavlovski et al.l [20061 ; ISchmidt et al]|2009l ). Consisten tly lower levels 



of \ < 0.15 were found in (adiabat ic) simulations at M s ~ 1 (jPouquet et al 



119911 : iPorter et al.lll994L Il999l . l2002h . As the strength of magnetic field fluctu- 



ations climbs up to equipartition with turbulent kinetic energy in our sequence 
of PPML simulations with Ma = 10, 3, and 1, the average level of x drops 
from 0.5 to below 0.3 for the sonic Mach number fixed at M s = 10 (see also 
Boldvrev et"aDl2002al lbh. 



To illustrate the effects of dynamic alignment in magnetized supersonic 
flows, in Fig. lc we show the probability density functions (pdfs) of the cosine 
of the alignment angle, cos# = u • ~B/Vu 2 B 2 , for three 512 3 PPML simulations 
at M s = 10 and Ma = 10, 3, and 1. In the most super-Alfvenic case at Ma = 10, 
the alignment is rather weak. It gets substantially stronger at Ma = 1, when 
the equipartition of turbulent magnetic and kinetic energies is reached. Fig, Id 
shows that the pdf of the alignment angle is well converged already at resolution 
of 512 3 in a series of representative numerical experiments with Ma = 3 and 
grid resolutions of 256 3 , 512 3 , and 1024 3 . 

Let us get back to Fig. lb and look at the shape of the compressional-to- 
solenoidal ratio m more detail. Since the two non-magnetized (PPM) simu- 
lations differ only in Xf and in the grid resolution, time-average spectra in Fig. lb 
indeed capture the effects of large-scale forcing that operates at kj /k m i n E [1, 2] 
in all cases shown. If one assumes that a small negative slope x{k) ~ A: -01 
seen at k/k m i n £ [10,60] (which can be traced to some extent in all simula- 
tions presented in this figure, see also Fig. le) is real, i.e. forms as a result of 
self-organization in supersonic turbulence, then it would seem that an isotropic 
forcing with Xf £ [0-6, 0.7] would be an optimal choice for M s ~ 6. At the same 



3 http: //lca.ucsd. edu/projects/enzo 
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time, the green line in Fig. lb shows that the effects of enforcing Xf = at 
kf j k m i n G [1,2] are felt at least up to k = I6k min , which is substantially larger 
than kf. A similar decline towards smaller wavenumbers is seen in the blue and 
pink curves, corresponding to Ma = 10 and 3, respectively. The trans- Alfvenic 
run at Ma = 1 shown in black does not have that feature. 

To better recognize the slope in x(k) ~ k~° , we replot the two 1024 3 results 
in a log-log plot, where it can be seen better (Fig. le). Is this slope real? Does it 
depend on the forcing? Is it related to specifics of numerical dissipation at small 
scales? To address these questions, we explored results of an MHD simulation 
of decaying turbulence, where the effects of continuous driving are minimized. 
This 1024 3 PPML simulation was carried out as part of KITP07 code comparison 
project^ The simulation follows a free decay of turbulence from a developed 
statistical steady state with M s ~ 10 and Ma ~ 10 down to M s ~ 2. Figure If 
shows x{k) for 10 flow snapshots equally spaced in time. As can be seen, with no 
forcing, the slope of —0.13 ± 0.03 is clearly present even though in this case we 
only have instantaneous power spectra and the data are rather noisy. In Fig. lg 
we show x(k) obtained for the 5th snapshot at 512 3 and 1024 3 to illustrate grid 
convergence with PPML. We also compared PPML results with those from three 
other popular numerical methods an d implementations fo r compressible ideal 
MHD (ZEUS, FLASH, RAMSES, see lKritsuk et alJ^OOQd ) at a grid resolution 



of 512 3 zones and found an excellent agreement. It seems that the slope is indeed 
real and does not depend on forcing or numerical dissipation. We derived a 
simple fitting formula for x(k) in this decaying turbulence model, assuming a 
fixed slope of —0.13 and fitting a normalization constant x(k/k m i n = 1) against 
Mg , as shown in Fig. lh. The approximation valid in the inertial range for rms 
turbulent Mach numbers M s £ [1, 10], 

X (k) « (0.44 + 0.004M, 2 ) {k/k min )- QA \ (1) 

is not unique, but still can be used as a guide in future experiments with opti- 
mized forcing for supersonic super-Alfvenic turbulence. 



3. Conclusion 

We explored the problem of forcing optimization to maximize the effective 
Reynolds numbers achieved in simulations of supersonic molecular cloud tur- 
bulence. Our results demonstrate that for low sonic Mach numbers and/or 
Alfvenic Mach numbers Ma ^$ 1 solenoidal forcing is reasonable. At the same 
time, for sonic and Alfvenic Mach numbers above 3, solenoidal forcing is inap- 
propriate, a proper mixture of dilatational and solenoidal modes is required for 
optimal representation of the inertial range. A purely compressional driving is 
impractical as it does not allow to study the scaling properties of fully developed 
supersonic turbulence. 
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